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A general relativistic version of the Euler equation for perfect fluid hydrodynamics is applied to 
a system of two neutron stars orbiting each other. In the quasi-equilibrium phase of the evolution 
of this system, a first integral of motion can be derived for certain velocity fields of the neutron star 
fluid including the (academic) case of co-rotation with respect to the orbital motion (synchronized 
binaries) and the realistic case of counter-rotation with respect to the orbital motion. The velocity 
field leading to this latter configuration can be computed by solving three-dimensional vector and 
f*** ' scalar Poisson equations. 

G\ '. 

Q\ . PACS number(s): 04.25.Dm, 04.30.Db, 04.40.Dg, 95.30.Sf, 97.60.Jd 

T— I ' 
^— > . 

o . 

O ; I. INTRODUCTION 

: 

Considerable efforts by many groups in the world are currently devoted to the computation of the gravitational 
y—i . radiation from binary neutron star coalescences (see e.g. jjj or [|| for recent reviews). These phenomena constitute 
one of the most promising sources of gravitational radiation for the interferometric detectors GEO600, LIGO and 
VIRGO currently under construction ||,[| . Basically two different approaches are used to tackle this problem: 

Al high-order Post-Newtonian analytical calculations in the point-mass limit for the two coalescing objects (see |^] 
for a review); 

t-h ; ; 

A2 hydrodynamical numerical simulations which treat the neutron stars as perfect fluid balls. In this latter category, 
different methods, based on different approximations, can be distinguished: 

O ! . rr- 

A2a Newtonian afflne approximations which consist in modeling the stars by triaxial ellipsoids, thereby reducing 
the dynamical degrees of freedom of a star to a finite number and leading to ordinary differential equations 
for the evolution, instead of partial differential equations f^-|| ; 

A2b Post-Newtonian afflne approximations, at the 1-PN order |9|-|Tl| or at the 2-PN order (l^l ; 

A2c Newtonian hydrodynamical simulations, either based on finite difference methods |l^-[l7j or on the Smooth 
Particles Hydrodynamical (SPH) method p8|-pl| (see ref. |l7j] for a comparison of various codes); 

A2d Post-Newtonian (at the order 1-PN) hydrodynamical simulations with a finite difference method @,^3| ; 

A2e fully relativistic hydrodynamical simulations within the 3+1 formalism of general relativity and using the 
approximation of a conformally flat spatial 3-metric along with no gravitational field dynamics p4|-p8t (see 
also H|). 

The analytical Post-Newtonian approach (Al) allows to compute the evolution of the binary system from an 
arbitrary early stage, when the separation between the two components is large, up to the rapid inspiral phase driven 
by the rapid loss of orbital energy by gravitational radiation. This approach breaks down when finite size effects 
(tidal forces, disruption of one of the stars) become important, i.e. during the coalescence phase. This final stage can 
be studied only by means of the numerical hydrodynamical methods (A2). But in this case one faces the problem of 
the initial conditions. Indeed, due to the limitation of computer resources, the initial conditions cannot be set when 
the separation between the two stars is much larger than their radii: this would require a prohibitive number of time 
steps for the evolution codes: the time to coalescence increases with the fourth power of the initial separation ao 
between the two objects. In practice, all the fully hydrodynamical computations listed in (A2) have been performed 
with ao set to at most a few times the stellar radius R: cio ~ 5R in ref. ao ~ 4i? in refs. [Q, |1S|| , ao ~ 3R in 
refs. |14|, |lq|, |i~7| ], ao = 2R (!) in ref. p2| . For the calculations employing the afflne approximation instead of the 



full hydrodynamics, the initial separation is taken to be somewhat larger: ao ~ 5R in ref. |9|, ao — 15 R in ref. [ p.2| 
Now at such small separations, two effects are important: (i) tidal forces, i.e. the influence of the gravitational field 
of star 1 (resp. 2) on the internal structure of star 2 (resp. 1), and (ii) general relativity. 
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Tidal effects have not been taken into account in the computation of the initial conditions of most of the fully 
hydrodynamical studies listed in (A2), the only exceptions being the works by Nakamura & Oohara jL3|, Rasio & 



Shapiro ]2Cj ] and Baumgarte et al. [|26|-28|, all in the case of synchronized binaries, i.e. when the stars have zero spin 
in the frame co-moving with the two centres of mass. However this rotation state is unrealistic. Indeed it can be 
seen ||[50) that the neutron star matter shear viscosity is too small and the binary evolution too rapid to lead to a 
synchronization of the spin periods with the orbital period^]. In other words, the inspiral of the binary system can 
be described in terms of perfect fluids and, in first approximation, all the forces acting on the fluid are gradients of 
some scalar potentials (gravitational force, gravitational radiation reaction force)^, so that Kelvin's theorem applies: 
the circulation of the fluid velocity (with respect to some inertial frame) on any closed contour comoving with the 
star (e.g. the stellar equator) is conserved. For each star, the circulation around the equator is roughly C ~ 2Auj 
where A is the area of the equatorial cross section of the star and u> its mean angular velocity with respect to some 
inertial frame 1Z: lu — f2 rb + ^spin where tt OI b is the angular orbital velocity of the system with respect to 1Z and 
f^spin is the rotation angular velocity of the star with respect to the co-orbiting frame. Since the variation of A is 
small during the evolution to the coalescence, the conservation of C is equivalent to the conservation of u>. When the 
separation is large f2 rb is negligible and u> is equal to the rotation rate of the star. When the separation is that of 
the initial conditions of the numerical computations (A2) (a ~ 50 km) f2 or b ~ 2x 10 3 rads _1 , which is much larger 
than lu, except for neutron stars rotating initially at millisecond periods. Consequently, one must have 

^spin — ^orb ■ (1) 

in order to have a constant circulation. We call configurations obeying ([j]) counter-rotating configurations. They 
represent realistic initial conditions for neutron star binary coalescence. 

Some of the fully hydrodynamical calculations listed in (A2) employ ([j]) as initial conditions . But 

none of them take into account the tidal effects: the stars are taken to be either spherical p5|,pl[ or axisymmetric 
1 14 1 (as mentioned above, the only computations with self-consistent initial conditions concerns synchronized binaries 
[13,|0 26-2f|, for which f2 SP j n = 0). To date, the only self-consistent initial conditions obeying (|]) have been computed 
by Bonazzola & Marck |32|. As can be seen in Fig. 3 of ref. the tidal deformation is quite important when the 
separation is ao ~ 3 R. However, these initial conditions have not been employed in evolution calculations. 

It must be noticed that some of the studies performed in the affine approximation (A2a) make use of self-consistent 
counter-rotating initial conditions. They correspond to irrotational Darwin- Riemann ellipsoids ||]8| or (in the ap- 
proximation of a large separation) irrotational Roche- Riemann ellipsoids gj. 

As regards general relativistic effects, the often used Newtonian approximation [items (A2a) and (A2c)] is very 
crude, in particular for the neutron star internal structure: let us recall that for a typical 1.4 Mq neutron star, 
the central value of the metric coefficient goo is around 0.4, which shows that even a (first order) Post-Newtonian 
approximation is not sufficient for describing these objects. 

The purpose of the present article is to give a method for computing self-consistent (i.e. including the tidal and 
rotational distortion) and realistic (i.e. obeying ([!])) initial conditions for binary neutron stars in the framework of 
the full general relativity. Therefore, this work can be conceived as the extension to general relativity of Bonazzola & 
Marck results |§]. 

The envisaged problem can be decomposed in two parts: (i) the computation of the gravitational field (i.e. the 
spacetime metric) generated by the two stars and their motion and (ii) the computation of the stellar structure 
(density distribution, velocity field,...) in that gravitational field. Part (i) is the main topic of numerical relativity and 
can be, at least in principle, be achieved by means of the classical 3+1 formalism (see e.g. j}3| and |Q). This paper 
focuses on the determination of the stellar structure. For this purpose, we consider that in the vicinity of the searched 
initial conditions, the system evolves along a sequence of equilibrium states. We obtain a first integral of motion for 
certain classes of velocity field inside the neutron stars, including the co-rotating and the (realistic) counter-rotating 
cases. 

The plan of the paper is as follows. Sect. || translates the basic assumption of quasi-equilibrium in geometrical 
terms (a spacetime symmetry) which leads to the definition of a privileged observer (the "co-orbiting" observer). The 



(relativistic) Euler equation for the fluid velocity is then derived in the frame of that observer (Sect. III). Necessary 



1 As shown by Kochanek j5j, this conclusion remains valid even if one takes into account the much higher effective viscosity 
arising when the neutron star's solid crust enters the plastic flow regime. 

2 this is not true for the so-called "gravitomagnetic" force; this latter induces some circulation of the fluid, as studied recently 
by Shapiro However this effect is important only for neutron star - black hole binaries, with a maximally rotating black 
hole. 
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and sufficient conditions to get a first integral of the Euler equation are given in Sect. IV. This first integral is trivially 
obtained in the (unrealistic) case of co-rotating stars. The astrophysically relevant case of counter-rotating stars is 
presented in Sect. M. A method of resolution is discussed in Sect. VI. 



II. SPACETIME SYMMETRY AND CHOICE OF COORDINATES 



A. Quasi-equilibrium hypothesis 

When the separation between the centres of the two neutron stars is about 50 km (in harmonic coordinates) the time 
variation of the orbital period, P or b, computed at the 2-PN order by means of the formulas established by Blanchet et 
al. Q is about 2%. The evolution at this stage can thus be still considered as a sequence of equilibrium configurations. 
Moreover the orbits are expected to be circular (vanishing eccentricity) , as a consequence of the gravitational radiation 
reaction |36| . In terms of the spacetime geometry, we translate these assumptions by demanding that there exists a 
Killing vector field l a which is expressible as 

l a = k a + Q olh m a , (2) 

where f2 or b is a constant, to be identified with the orbital angular velocity with respect to a distant inertial observer, 
and k a and m a are two vector fields with the following properties. k a is timelike at least far from the binary 
and is normalized so that far from the star it coincides with the 4-velocity of inertial observers. m a is spacelike, 
has closed orbits, is zero on a two dimensional timclike surface, called the rotation axis, and is normalized so that 
V Al (m p m p )V At (m cr m tT )/(4TO^TO 1/ ) tends to 1 on the rotation axis [this latter condition ensures that the parameter 
4> associated with m a along its trajectories by m a = (d/d(f)) a has the standard 2ir periodicity]. Let us call l a the 
helicoidal Killing vector. We assume that l a is a symmetry generator not only for the spacetime metric gap but also 
for all the matter fields. In particular, l a is tangent to the world tubes representing the surface of each star, hence its 
qualification of helicoidal (cf. Figure yj). 

The approximation suggested above amounts to neglect outgoing gravitational radiation. For non-axisymmetric 
systems — as binaries are — it is well known that imposing l a as an exact Killing vector leads to a spacetime which 
is not asymptotically flat |3^ ]. Thus, in solving for the gravitational field equations, a certain approximation has to 
be devised in order to avoid the divergence of some metric coefficients at infinity. For instance such an approximation 
could be the Wilson & Mathews scheme |38| that amounts to solve only for the Hamiltonian and momentum constraint 
equations. This approximation has been used in all the fully relativistic studies to date p4]-|28| and is consistent with 
the existence of the helicoidal Killing vector field (Q) . Note also that since the gravitational radiation reaction shows 
up only at the 2.5-PN order, the helicoidal symmetry is exact up to the 2-PN order. 



B. 3+1 Foliation of spacetime 



For the considered problem, two types of coordinates can be envisaged: "non-rotating" coordinates (t,x z ) which 
are Minkowskian at infinity, so that k a is the first vector of the natural basis corresponding to these coordinates and 
"co-rotating" coordinates (t',x z ) so that l a is the first vector of their natural basis. There is a lot of ways to do this. 
We choose both coordinate systems so that the hypersurfaces t = const and t' = const coincides and are maximal 
spacelike hypersurfaces. More precisely, we suppose that there exists a slicing of spacetime by a family of spacelike 
hypersurfaces (S t ) so that (i) each E t is spacelike and (ii) m a is tangent to £(. 



1. Non-rotating coordinates 

On each S t , we choose a system of Cartesian coordinates x % — (x, y, z), such that (t, x l ) is a system of spacetime 
coordinates satisfying to 

*■-(!)' _ _ < 3 » 

~»(s) " + "(£)"■ (4) 
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The lapse function N and shift vector N a associated with these coordinates are defined by 

k a = Nn a - N a and n^N" = , 
where n a is the future directed unit 4- vector normal to the hypersurface E t . 



(5) 



2. Rotating coordinates 



We call rotating coordinates any coordinate system x % on S t such that (t! — t, x l ) is a spacetime coordinate system 
satisfying to 




In other words, the lines x l = const are the trajectories of l a . This latter being a Killing vector, this means that 
t' is an ignorable coordinate for such systems. In numerical studies, we will use these types of coordinates to reduce 
the a priori 4-D problem to a 3-D one. In practice, three rotating coordinate systems can be used: one centered on 
each star, to describe properly the hydrodynamics, and a third one centered on the rotation axis, to describe the 
gravitational field. The lapse function and shift vector associated with rotating coordinates are immediately deduced 
from Eqs. (|) and (g) which result in 

l a = Nn a - B a , (7) 

with 



B a := N a - 0, rb m a . (8) 

Since m a is parallel to S t , B a is indeed the shift associated with rotating coordinates (cf. Figure |l|). Note that 
rotating and non-rotating coordinates have the same lapse N for they define the same spacetime foliation. 

The Killing equation V a lp + V^/q = 0, once projected onto Ej, leads to the following relation between the extrinsic 
curvature tensor K a p of the hypersurfaces St and the derivatives of the shift vector B a : 

2NK afj = -V a B p - V f3 B a - 2n a nf3n^V^N (9) 

or, equivalently, 

2NK a0 = -V a Bf) - VpB a , (10) 

where V a stands for the covariant derivative associated with the 3-metric h a p induced by g a p onto the hypersurfaces 
St. Note that K a p is linked to the covariant derivative of n a by the formula 

V l3 n a = -K a0 -V a (\nN)n l 3 . (11) 
In the following an extensive use is made of this relation, without explicitly mention it. 



C. The co-orbiting observer 



Let us call the co-orbiting observer the observer O whose world lines coincide with the trajectories of the symmetry 
group when these latter are timelike, which encompasses the region occupied by the two stars. O's 4- velocity v a can 
be written: 

u Q = e-*r, (12) 

where $ is a scalar field which is uniquely specified by the normalization relation v^v^ = — 1. It can be seen easily 
that <!> is related to the lapse N, the shift vectors B a and N a , the azimuthal vector m a and the orbital angular 
velocity fl OY h by 

e 2* = iV 2 - B„B^ = N 2 - n 2 olh m^m» + 2Q OTh m„N^ - . (13) 
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Let q a p be the projector onto the 3-planes II orthogonal to v a : 

q a /3 ■= g a /3 + v a vp . (14) 

The kinematics of the observer O is entirely specified by the Ehlers decomposition |3{| of the covariant (with respect 
to gap) derivative of v a : 



where 



is the rotation 2- form of O, 



is the expansion tensor of O and 



V(3V a = Lu al3 + 9 a/3 - a a Vj3 , (15) 

w Q /3 := q^q^V^ (16) 

9 aP := q a "q u V (v w m) (17) 

a a := v^V^a (18) 



is the 4~ acceleration of O. The property (|l2|), namely that n" is collinear to a Killing vector, means that v a is an 
isometric flow |39f| and leads to 

Q /3 = (19) 

and 

a a = V a $ . (20) 

Equation ( |l9| ) shows that i> Q is a rigid flow. 

The three-dimensional vector space II represents the local rest frame of O. Note that since O is rotating (ui a p ^ 0, 
see below), II is not integrable into global 3-surfaces. q a p is the (positive definite) metric tensor induced by g a p on 
II. We can introduce the alternating tensor within II as 

Cq/37 := V^Cftapj , (21) 

where e a f3-ys is the spacetime alternating tensor associated with the spacetime metric g a p- The rotation 2-form of O 
is fully specified by its dual within II: u a p — —uj^t^ap, where 



oj := —-e^u 



2 - -mu 2 



e a ^V„v v . (22) 



Note that the Raychaudhuri identity for the flow v a reduces to a simple relation between the norm of tu a and the 
Laplacian of 

w M u" = -V M V$ + X -R^v v , (23) 

where R a p is the Ricci tensor of the metric g a p. 

By means of Eqs. (|l2|), (f7|), (||) and the Killing identity V a lp + V pl a = 0, the rotation vector ( p2| ) can be expressed 

as 

w° = y e a ^ [O or bV^m„ - + 2(fi orb m p - jV M )V v IniV] . (24) 
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III. RELATIVISTIC EULER EQUATION IN THE ROTATING FRAME 



A. Fluid motion 

As stated in the introduction, the matter constituting the neutron stars can be considered as a perfect fluid, so 
that its stress-energy tensor writes 

T al3 = {e+p)u a u p +pg ap . (25) 
The fluid 4-velocity u a can be decomposed orthogonally with respect to the rotating observer O as follows: 

u a =T(V a + v a ), (26) 

where T is the Lorentz factor 

r := , (27) 

and V a is the fluid 3-velocity with respect to O: 

V a := . (28) 

V a belongs to II and is the fluid velocity as measured by the observer O (i.e. with respect to Cs proper time). As 
an immediate consequence of u^u^ = — 1, one has the usual relation between T and V a : 

r = (i-T^v^r 1/2 . (29) 

To a very good approximation the (cold) neutron star matter equation of state (hereafter EOS) can be considered 
as barotropic: e = e(n) and p — pin) where n is the proper baryon density. It is then worth to introduce the logarithm 
of the ratio of enthalpy per baryon and the baryon mean rest mass tub by 

H:=ln( e -±1) , (30) 
\m B n J 

which we call the log- enthalpy, to re-express (V Q p)/(e + p) as a gradient of a scalar: indeed, by virtue of the First 
Law of Thermodynamics, the following identity holds for any barotropic EOS: 

^ - V aJ ff . (31) 

e + p 

Then the fundamental energy-momentum conservation equation 

V^r Q = (32) 

can easily be shown to be equivalent to the system of the following two equations 

V^™/ 1 ) = (33) 
u^ 1 V m m q + V a H + u^V^H u a = . (34) 



Note that Eq. (34) is nothing but the uniformly canonical equation of motion for a single-constituent perfect fluid as 
given by CarterQ. 

B. Baryon number conservation 

Inserting Eq. ( p6| ) into the baryon number conservation equation ( |33] ) and using the fact that v a is divergence-free 
[Eq. (|l9|)] leads to 

V M y /i + V , V A ,In(nr) = 0. (35) 



G 



C. Momentum conservation 



By projecting Eq. (Q) onto v a (i.e. along the world lines of the co-orbiting observer O), one obtains the relation 

V V p (H + $ + In T) = , (36) 

which can be considered as a relativistic generalization of the classical Bernouilli theorem, for it means that the 
quantity H + $ + lnT is constant along the fluid lines. 

By projecting Eq. (|3~i| ) perpendicularly to v a (i.e. onto the local rest frame of the co-orbiting observer 0), one 
obtains the relativistic version of the Euler equation for the fluid velocity with respect to O : 

FV^r + 2e a llu uj^V" + V Q $ - V^V^ V a + T- 2 \7 a H = . (37) 

At the Newtonian limit, the V^V ^V a gives the classical term (V-V)V and 2e a fll/ cu fJ -V l/ gives the Coriolis term 2u)xV, 
induced by the rotation of the observer O with respect to some inertial frame. The term V a H gives the classical 
pressure term. Finally reduces to the sum of the gravitational and centrifugal potentials [cf. Eq. (filf)] 

$ N ^ $ grav _I(fUxr) 2 , (38) 

^grav being defined so that the gravitational field writes g = — V<&g rav . 

In order to exhibit from Eq. ( p7|) a first integral of motion, we shall write as much terms as possible under the form 
of gradients. First, it can be seen easily that, similarly to the usual flat space formula, the following identity holds 

V»V^V a = e a ^(\7 AVrV" + V a f-V^vJ , (39) 

where we have introduced the curl of V a within the 3-space II: 

(V A V) a := e a M!/ VV . (40) 
Putting Eq. ( ^9|) into Eq. ( |37j ) and performing slight rearrangements results in 

v a (H + $ + inr) + r 2 {i ailv [(v a vy + 2cy<] v v + ^r?/v„$} = o , (4i) 



where 



is the projector onto the 2-space orthogonal to V a , i.e. orthogonal to the fluid lines with respect to O. Note that 
in the case where the fluid is at rest with respect to O, Pj 3 is not defined; however, the product V il V ix P ( J i which 
appears in Eq. ( f4l| ) remains well defined and is equal to zero. In the derivation of Eq. (|4l|), use has been made of 
Eq. (H) to replace the term V a (V^^/2) coming from Eq. © by r- 2 V Q lnr. 



D. Number of independent components 

From the fundamental equation V M T Ma = 0, which has a priori four independent components, we have derived 
two scalar equations [Eqs. (|35| ) and Eqs. (|3^)] and one vectorial equation [Eq. Jii])]. Equations ([36]) and ( [4l| ) are not 



independent: the former is a direct consequence of the latter, as seen easily by projecting (|4lj) along V a . Moreover, 
from its construction, pi] ) has only three independent components for it lies into the 3-planes II orthogonal to v a . 

We will take the scalar equation ( |35| ) and the vectorial equation lying in II ([y]) as the fundamental equations to 
be satisfied for our problem. 
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IV. CONSTRAINT ON THE VELOCITY FIELD AND FIRST INTEGRAL OF MOTION 



The only assumption underlying Eq. ( f4l| ) is that the observer O, with respect to which the fluid velocity V a is 
defined, has world lines parallel to the helicoidal Killing vector l a . Equation ([hy) is equivalent to the system 

e a »u [(V A VY + 2a/ 1 ] V + V^PJVv® = T- 2 V Q G (43) 

H + $ + lnT + G = const. , (44) 

where G is a scalar field defined at least in the stars' world tubes. The relation (Q) constitutes a first integral of 
motion of the system. 

A. The co-rotating case 

Equation (^) is trivially satisfied in the case V a — by taking G = const. The first integral of motion reduces 
then to 

H + $ = const. (45) 

This case is the co-rotating one (f2 S pm = 0) mentioned in the introduction; it corresponds to synchronized binaries, 
which are not expected to represent realistic close neutron star binaries, as discussed in Sect. |j. 

The first integral (|45|) follows directly from the fact that in the co-rotating case the fluid 4-velocity is parallel to a 
Killing vector pi]: V a = is indeed equivalent to u a = v a — e~®l a [cf. Eq. (26)]. The integral (45) is well known in 



the case of a single rigidly rotating star (see e.g. ]]42j). For the problem of the initial conditions of a binary coalescence, 
it has been used by Nakamura & Oohara (at the Newtonian approximation) and Baumgarte et al. [ p6[ - 

B. Formulation of the problem in the general case 

From now on, we suppose that V a ^ 0. By performing the vector product (with respect to tap-y) of Eq. ( f43| ) by 
V a , one can see easily that Eq. (|i^) is equivalent to the system 

^ V M G = (46) 
(V AV) a = -2oj a -e a ^V^^+{T 2 V ry V' 7 y 1 e a ^V^V v G + FV a , (47) 



where 



1 



F '■= T7T7Z V v [( V A V T + 2w 1 ■ (48) 



The gravitational field being given, the problem of getting a solution amounts to finding a vector field V a and a 
scalar field G such that the equations (|35|), (^) and (|47j) are satisfied. In Eq. (|35|), the scalar field n is that related 
to the gravitational field, V and G by the first integral ( El ) via the EOS n = n(H). More precisely, let us consider 
an iterative method for solving this problem. Let us suppose that at a given step, the gravitational field equations 
have been solved; the potential $ and the rotation vector Lu a are then known. The enthalpy H can be then deduced 
from the first integral (^Z|) , by taking for T and G the values at the previous step or making some extrapolation from 
a few previous steps. The baryon density n is computed from H by means of the EOS. The system of equations (|3^), 
(|h|) and (^) is then to be solved in V a and G. It is however not obvious that a solution exists in the general case. 
What can be said is that in the Newtonian and incompressible case, solutions do exist and are constituted by S-type 
Darwin- Riemann ellipsoids |43fl . 

V. THE COUNTER-ROTATING CASE 
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A. Definition 



Let us focus on the interesting case of counter-rotating binaries. The concept of counter-rotation has to be defined 
in the relativistic framework. We shall define it by requiring V a + 1 and the scalar field G introduced in Eqs. (|43|)-(|44|) 
to be constant: 

G = const. (49) 

This definition is motivated by the fact that at the Newtonian limitj^] it implies 1/2 (V A V) a = — u) a , which is the 
definition (|l|) of counter-rotation [cf Eq. (p4])j. 

With the choice (|49|), Eq. ( |4^ ) is trivially satisfied and Eq. ( [47|) becomes 

(V A V) a = -2uj a - e^V^Vv® + FV a . (50) 



B. 3+1 decomposition 

From the numerical point of view, it is desirable to reduce the problem to the resolution of three-dimensional 
equations. Now Eq. (|50| ) involves 4- vectors: even if V a is spacelike and belongs to 3-planes orthogonal to v a , due 
to the rotation of this latter, there exists no coordinate system in which V a would have only three non-vanishing 
comp onen ts. Therefore, we choose to recast Eq. ( |50| ) according to the 3+1 foliation of spacetime introduced in 



Sect. [IB. In this manner, we will consider only 3- vectors belonging to the spacelike hypersurfaces £(. The first step 
is to introduce the orthogonal decomposition with respect to S t of the fluid velocity V a with respect to the co-orbiting 
observer: 

V a = W a + Zn a , (51) 

where 

W a = /i Q ^ and Z = -n M V" , (52) 

where h a p := g a p + n a np is the orthogonal projector onto E t , or equivalently, the 3-metric induced by g a g in E t . 
Due to the orthogonality relation v^V^ = 0, the scalar Z is not independent from W a : by inserting Eqs. (12) and (Q) 
into v^V^ — 0, one gets 



Z = ~B IA W> l = ~B i W i (•+!! 



In the last part of this equation, we have introduced Latin indices, which range from 1 to 3, whereas the Greek indices 
range from to 3. We will systematically do this in the following for all the tensor fields that lie in E t , such the 
vectors B a and W a . In this way the three-dimensional character of the equations will clearly appear. 

The curl on the left-hand side of Eq. (^fj| ) is defined within respect to the alternating tensor e a ^ 7 , which is neither 
parallel nor orthogonal to the (T, t ) foliation. Let us introduce instead the alternating tensor e Q/37 within the space 
(E f , hap) by 

e^T : = n^Pi . (54) 

Inserting the identity e a / 3 T 5 = —^ n [ a ^0i s ] into the definition ( |2l"| ) of e a ^T, we arrive at the following expression of e^ 1 
in term of l a ^\ 

= e -*(jv^ - n a B^i - n^B^ 01 - n?B^ al3 ) . (55) 

Besides, we can express the four-dimensional covariant derivative of V a which appears in the curl of V a in terms 
of the three-dimensional covariant derivative of W a with respect to h a p ■ 



Details about the Newtonian limit will be presented in Sect. VD 
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VaVJj = V a Wp ~ n a n^ M Wp - K^Wnp + n V a Z - ZK a p - Zn a V > In AT 



(56) 



A useful formula is that which gives the derivative along n a of any tensor field X a (i) which lies in Ej and (ii) 
which respects the heloicoidal symmetry (i.e. C\X a = 0): 



This formula can be used to express the derivative of W a along n a which appear in the right-hand side of Eq. (]5 
We can also apply it to B a and get 



(57) 



n"Vf,B a = -K^B" 



B^ ~ 



(58) 



By combining Eqs. 



(VAV) 



a g — ^ ^ a t lu 



j); ( p5| ) an d (|57j), we arrive at the 3+1 decomposition of (V A V) a : 

B„ 



N 



AV„ir„ + (B°V a W v - W a V a B„) - Z^V„N + B v V, t Z - 2K°W*B V 



N 



_ e -9 n * e^BxVaWu 



>„wu- (59) 

The term with ( a ^ v in factor is parallel to the hypersurface E t (because i afJ,L/ is), whereas the second term is along 
n a . 

Similarly Eqs. (§|), (H), (@), © and (||) leads to the 3+1 splitting of the rotation vector u a : 

N . 



e -2$ gx\iv 



-y^B v - £ M V„Ar + K a B a B v 



-2<I> 



i a e x ^B x V^B v . 



(60) 



We also need to perform the 3+1 decomposition of the second term on the right-hand side of Eq. (|50|). The result 



is 



(61) 



Thanks to Eqs. (|59j), ( pPy and (|6l|), the orthogonal projection of Eq. ( |50| ) onto the hypersurfaces E t is straightforward 
and leads to the three-dimensional equation 



NVjW k + ^ (B l ViW k - W l ViB k ^ - Z^V k N + B k VjZ - 2K j l W i B k 



JVVjB fc + 2B 3 W k N - 2K j 1 B t B k ) - AW^V^ - ZBf7 k $ + — WjBkB'Vj* 



(62) 



The baryon number conservation equation (p5|), once recast in terms of three-dimensional quantities with the help 
of Eq. J56|)) writes 



ViW % + W' l V l ln(ATra) + —B i V l \n{ZTn) = . 



(63) 



A boundary condition on W l can be derived by multiplying this equation by n and setting n = (definition of the 
star's surface) into the result. One obtains, using Eq. (p3|), 



N 2 



V» n 







(64) 



surface 



The equations to be solved are the three-dimensional vector equation (|62j) and the scalar equation (63), altogether 
with the boundary condition (|64|). Note that once Eq. (^) is satisfied, the other part of the four-dimensional 
equation ( j50|), n amely the part along n a , is automatically fulfilled. Indeed, the projection of Eq. ( |50| ) along n a leads 
to [cf. EqsJ^, ©, © and (§§] 

e^V^ = e-*eV k BitjB k - e y7c £^ V fe $ + F , (65) 
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which is nothing else than the orthog onal p rojection of Eq. ( |62| ) onto B % , 

Referring to the discussion in Sect. Ill D , we conclude that if a 3- vector W l of St obeying Eqs. ( p2|) and (|63] ) can be 
found, the problem is completely resolved. Note that the baryon density n which appears in Eq. (|6^) is given via the 
EOS by the log-enthalpy H, itself being fully determined (for a fixed gravitational field) by W l via the first integral 
of motion (H) , which becomes in the counter-rotating case 



H + $ + lnr = const. , 



where [cf. Eqs. (|29|) and (|T 



r = (i - w,w l + z 



2N-1/2 



(66) 



(67) 



C. Formulation in terms of Poisson equations 



The equations to be solved, namely Eqs. ( |62| ) and (|63|), can be recast into Poisson equations by looking for solutions 
under the form 



(68) 



where ip is a scalar field and A 1 is a 3-vector of S t , which without any loss of generality can be taken to be divergence- 
free : 



ViA 1 = . 

This latter property implies that the curl of W l is related to the Laplacian of A % by 

e ijk VjW k = -V j V j A i +R l ,A j , 



(69) 



(70) 



where Rij is the Ricci tensor of the 3-metric hij of S t . Inserting Eq. ( |70[ ) into Eq. (62) leads to the following vector 
Poisson equation for A 1 : 



' l\/ Ihlh.) ■ U ,V;«I- ' 



N 



N 



N J 



(B l ViW k - W l ViB k >) - Z^-VkN + B k ij 3 Z - 2K J l W t B k j 



N 



The divergence of W l evaluated from Eq. (pq) is 



ViW 1 = V,VV + h ijk R kHj A 
2 



W l + R i „A j . (71) 



(72) 



where Rkiij is the Riemann tensor associated with the 3-metric hij. By virtue of the symmetry properties of the 
Riemann tensor in three dimensions, the last term on the right-hand side of Eq. ([72]) vanishes identically, so that the 
divergence of W % is simply the Laplacian of tp and Eq. (p3f) becomes 



ViVY = -W*Vi ln(ATn) - jjB^i In(ZTn) . 



(73) 



D. Newtonian limit 



In the Newtonian limit, $ takes the form ( p8| ) and uo a becomes [cf. Eq. (|24|)1 

a Newt. in a\ 

0J ► ^orb ■ (74) 

The rotating-coordinate shift vector B a reduces to B a = —Vt or \>m a [cf Eq. (B)], so that Eq. (fn|) becomes 
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AA = 2 orb) (75) 

whose divergence-free solution is 

A= \(e z xrf J2 orb , (76) 



where e z :— fi rb/^orb- Finally, Eq. (73) becomes 

A^ = -V ■ Vlnn . (77) 

Once this equation is solved, the fluid velocity field with respect to the co-orbiting observer is computed by taking 
the curl of © : 

V = -Oorb x r + VV> . (78) 



This is the solution obtained by Bonazzola & Marck 32 



VI. DISCUSSION 
A. Iterative method of resolution 

The resolution of the problem amounts to solving the vector Poisson equation ( |7l| ) for A % and the scalar Poisson 
equation (|73"|) for ip, with the boundary condition ( pi]) at the surface of each star. These equations involve Laplacian 
with respect to the curved 3-metric hij, so that even if the right-hand of the equations is supposed to be known (e.g. 
from a previous step in an iterative method), the numerical solution is not straightforward. A technique which has 
shown to be successful consists in introducing on S t a flat 3-metric hij and decomposing the operators into flat-space 
ordinary Laplacians plus curvature terms Q,[l5|. 

With this technique, the following iterative method can be envisaged to get counter-rotating binary neutron star 
configurations. The starting point of the procedure can be very crude approximations such as constant density 
spherical stars with W % = and a flat spacetime metric. At a given step, the gravitational field equations are to be 
solvedR. leading to new values for the functions N, hij, Kij, B l and $ [via. Eq. (|l^)]. The first integral of motion 
Eq. ( pq ) yields then to the value of the log-enthalpy H throughout the stars. The Lorcntz factor T which appears in 
the first integral is to be evaluated by inserting the previous step values of W l in Eq. ( |67| ) . From H the baryon number 
density n is computed by means of the EOS and inserted in the right-hand side of the scalar Poisson equation J75|) . 
In this right-hand side, as well as in the right-hand side of the vector Poisson equation (|7l]), the value of W % is to be 
taken from the previous step. This is of course also the case of the functions of W l : Z [deduced from W l by Eq. J53|)] 
and F [deduced from W l by Eq. (|4S|)]. The Poisson equations (|7^) and ( |7l| ) are then to be solved respectively for ip 
and A 1 . The solutions are obtained up to the addition of harmonic functions. These latter are determined in order 
that the boundary condition ( |64] ) is fulfilled. The vector field W l is deduced from the obtained values of ip and A % 
via Eq. (|6^) and a new iteration may begin. 

Of course, we do not have any theorem about the convergence of this iterative procedure. All that we can say 



is that similar schemes have been applied successfully to the computation of axisymmetric J42j and triaxial [|41 45 
models of single neutron stars and that in the axisymmetric case, a rigorous proof of convergence has been recently 
given by Schaudt & Pfister 

B. Conclusion 

Before the inner most stable orbit is reached, the evolution of a binary system of neutron stars can be approximated 
by a sequence of quasi-equilibrium con figur ations. For each of these configurations, the spacetime possesses the 



helicoidal symmetry discussed in Sect. II A. The hydrodynamical part of the problem is then trivial in the case 



4 Although the gravitational field equations are not discussed in the present paper, let us note that thanks to Eq. ([U]), the 
momentum constraint equation can be expressed as a three-dimensional vector Poisson equation for B % . 
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of synchronized binaries, because of the existence of the first integral of motion (p5[), which means that once the 
gravitational field is known, the matter distribution in the stars is obtained immediately. However realistic neutron 
star binaries on the verge of coalescence are not synchronized but rather in counter-rotation. In this case, the velocity 
field inside the stars with respect to the co-orbiting observer is not zero and has to be computed so that the Euler 
equation (|37]) is satisfied. We have presented a formalism which reduces the problem of finding this velocity field 
to the resolution of three-dimensional scalar and vector Poisson equations. We are currently applying the numerical 
techniques we have recently developed for solving such equations with spherical-type coordinates J44| , [45"t in order to 
get numerical models. We will present the results of this work in a future article. 

The formulation presented in this article is independent of any prescription for solving the gravitational field 
equations (Einstein equations). It simply relies on the assumption of the spacetime helicoidal symmetry and can be 



used in conjunction with any set of equations for the gravitational field, such as the Wilson & Mathews' scheme |58 25 
or the 2-PN scheme recently proposed by Asada & Shibata |47|. Note that both schemes involve nothing else but 
the resolution of Poisson type equations, so that the method that we propose do not require a numerical technique 
specific to the hydrodynamical equations. 
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